#### "Economic Inequality. Income and political participation in authoritarian regimes" ####
# authors: "Pelke, Lars"
# date: 2019-06-17
# written under "R version 3.6.0 (2019-03-11)"
# Robustness test vote models 

# Structure #
## 1. Income Deciles 
## 2. All the Ginis measure
## 3. Civil Society variable without party-near organizations
## 4. Multiple Imputation of missings

#### Preliminaries ####

R.version$version.string

# clear workspace
rm(list=ls())

# set working directory

# loading packages

library(countrycode)
library(tidyverse)
library(viridis)
library(lme4)
library(sjstats)

###############################################################################################
###############################################################################################
#### 1. Income Deciles 
#### Supplementary Appendix E1 and E2


#### Import Data ####
wvsdata <- readRDS("data/wvsdata.rds") 

#### Reducing to authoritarian regimes ####

wvsdata <- wvsdata %>%
  filter(v2x_regime_amb<=4)

wvsdata <- wvsdata %>%
  mutate(country_year = stringr::str_c(country_name, year, sep=" "))

table(wvsdata$country_year, is.na(wvsdata$gini_mkt))

#### Building variables ####

#Replacing missings in Dependent Variables

wvsdata <- wvsdata %>%
  filter(A064>=0 | A066>=0 | A067>=0 | A068>=0 | A069>=0 | A070>=0 |
           A072>=0 | A073>=0 | A074>=0 |A075>=0 | A076>=0 | A077>=0
         | A078>=0 | A079>=0 | A099>=0 | A100>=0 | A101>=0 | A102>=0
         | A103>=0 | A104>=0 | A105>=0 | A106>=0 | A106B>=0 | A106C>=0 )%>%
  mutate(civilsociety_b = if_else(A064==1, 1,
                                  if_else(A066==1, 1,
                                          if_else(A067==1, 1, 
                                                  if_else(A068==1, 1,
                                                          if_else(A069==1, 1, 
                                                                  if_else(A070==1, 1,
                                                                          if_else(A072==1, 1, 
                                                                                  if_else(A073==1, 1,
                                                                                          if_else(A074==1, 1,
                                                                                                  if_else(A075==1, 1,
                                                                                                          if_else(A076==1, 1, 
                                                                                                                  if_else(A077==1, 1,
                                                                                                                          if_else(A078==1, 1,
                                                                                                                                  if_else(A079==1, 1, 
                                                                                                                                          if_else(A099>=1, 1,
                                                                                                                                                  if_else(A100>=1, 1, 
                                                                                                                                                          if_else(A101>=1, 1,
                                                                                                                                                                  if_else(A102>=1, 1,
                                                                                                                                                                          if_else(A103>=1, 1,
                                                                                                                                                                                  if_else(A104>=1, 1, 
                                                                                                                                                                                          if_else(A105>=1, 1,
                                                                                                                                                                                                  if_else(A106>=1, 1,
                                                                                                                                                                                                          if_else(A106B>=1, 1, 
                                                                                                                                                                                                                  if_else(A106C>=1, 1, 0)))))))))))))))))))))))))

table(wvsdata$civilsociety_b)

# Dependent variable without party organizsations 
wvsdata <- wvsdata %>%
  filter(A064>=0 | A066>=0 | A067>=0 | A069>=0 | A070>=0 |
           A072>=0 | A073>=0 | A074>=0 | A075>=0 | A076>=0 | A077>=0
         | A078>=0 | A079>=0 | A099>=0 | A100>=0 | A101>=0 
         | A103>=0 | A104>=0 | A105>=0 | A106>=0 | A106B>=0 | A106C>=0 ) %>%
  mutate(civilsociety_p = if_else(A064==1, 1,
                                  if_else(A066==1, 1,
                                          if_else(A067==1, 1, 
                                                  if_else(A069==1, 1, 
                                                          if_else(A070==1, 1,
                                                                  if_else(A072==1, 1, 
                                                                          if_else(A073==1, 1,
                                                                                  if_else(A074==1, 1,
                                                                                          if_else(A075==1, 1,
                                                                                                  if_else(A076==1, 1, 
                                                                                                          if_else(A077==1, 1,
                                                                                                                  if_else(A078==1, 1,
                                                                                                                          if_else(A079==1, 1, 
                                                                                                                                  if_else(A099>=1, 1,
                                                                                                                                          if_else(A100>=1, 1, 
                                                                                                                                                  if_else(A101>=1, 1,
                                                                                                                                                          if_else(A103>=1, 1,
                                                                                                                                                                  if_else(A104>=1, 1, 
                                                                                                                                                                          if_else(A105>=1, 1,
                                                                                                                                                                                  if_else(A106>=1, 1,
                                                                                                                                                                                          if_else(A106B>=1, 1, 
                                                                                                                                                                                                  if_else(A106C>=1, 1, 0)))))))))))))))))))))))

table(wvsdata$civilsociety_p)

number_observations <- wvsdata %>%
  count(country_year)


#Repalcing missings in Independent Variables

wvsdata <- wvsdata %>%
  filter(X001>=0)%>%
  mutate(sex= X001)

wvsdata <- wvsdata %>%
  filter(X003>=0)%>%
  mutate(age= X003)

wvsdata <- wvsdata %>%
  filter(X025>=0)%>%
  mutate(education= X025)

wvsdata <- wvsdata %>%
  filter(X028>=0)%>%
  mutate(employement_status= X028)

wvsdata <- wvsdata %>%
  filter(X047>=0)%>%
  mutate(income_deciles= X047)

wvsdata <- wvsdata %>%
  mutate(income_quintiles = ifelse(income_deciles==1, 1,
                                   ifelse(income_deciles==2, 1,
                                          ifelse(income_deciles==3, 2,
                                                 ifelse(income_deciles==4, 2,
                                                        ifelse(income_deciles==5, 3,
                                                               ifelse(income_deciles==6, 3,
                                                                      ifelse(income_deciles==7 , 4,
                                                                             ifelse(income_deciles==8, 4, 5)))))))))


wvsdata <- wvsdata %>%
  filter(X011>=0)%>%
  mutate(children= ifelse(X011>0, 1, 0))

wvsdata <- wvsdata %>%
  filter(X028>=0)%>%
  mutate(unemployed= ifelse(X028==7, 1, 0))

table(wvsdata$civilsociety_b)

# 35907 non member of civil society/ not belonged to civil society, 29471 member/ Belong to cso

# 40217 non member of civil society/ not belonged to civil society, 32846 member/ Belong to cso

number_observations2 <- wvsdata %>%
  count(country_year)

## Inspect which countries are delete because missing obserations at macro-level variables ##

setdiff(number_observations$country_year, number_observations2$country_year)

# Comment: [1] "China 1990"        "Croatia 1996"      "Georgia 1996"      "Hong Kong 2014"    "Jordan 2007"       "Mexico 1981"      
# [7] "Mexico 1990"       "Russia 1990"       "South Africa 1982" "South Korea 1982" are excluded from the dataset

###Data Manipulation of GINI measures, fill some missing value with last national observations

missing <- wvsdata %>% 
  group_by(country_year) %>% 
  summarise(sumNA = sum(is.na(gini_mkt)))
missing

wvsdata <- wvsdata %>%
  mutate(gini_mkt=ifelse(country_year=="Algeria 2014", 36.7, gini_mkt)) %>% # Algeria 2011
  mutate(gini_disp=ifelse(country_year=="Algeria 2014", 35, gini_disp)) %>%
  mutate(gini_mkt=ifelse(country_year=="Azerbaijan 2011", 38.3, gini_mkt)) %>% # Azerbaijan 2011
  mutate(gini_disp=ifelse(country_year=="Azerbaijan 2011", 30.3, gini_disp)) %>%
  mutate(gini_disp=ifelse(country_year=="Lebanon 2013", 38.5, gini_mkt)) %>% # Lebanon 2012
  mutate(gini_disp=ifelse(country_year=="Lebanon 2013", 36.5, gini_disp))

wvsdata <- wvsdata %>%
  mutate(gini_mkt=ifelse(country_year=="Haiti 2016", 57.6, gini_mkt)) %>%
  mutate(gini_disp=ifelse(country_year=="Haiti 2016", 53.7, gini_disp)) %>%
  mutate(gini_mkt=ifelse(country_year=="Nigeria 2012", 46.8, gini_mkt)) %>%
  mutate(gini_disp=ifelse(country_year=="Nigeria 2012", 44.1, gini_disp)) 

### delete Kuwait 2014, Uzbekistan 2011 and Libya 2014, no GINI information in a range of maximum 4 years before or after the country-year observations

wvsdata <- wvsdata %>%
  drop_na(gini_mkt)
rm(missing)

table(wvsdata$civilsociety_b)
# 36955 non member of civil society/ not belonged to civil society, 30677 member/ Belong to cso

#### Previous democratic experience ####

table(wvsdata$democratic_expericene) 
table(wvsdata$democratic_expericene_all)

democratic_exp <- wvsdata %>%
  filter(democratic_expericene_bi ==1) %>%
  group_by(country_name) %>%
  count()
# authoritarian regime that has previous democratic experience: Armenia (1990-1994), Belarus (1991-1996), Palestine (2004-2006), Russia (1992, 1993), 
# Thailand (1998-2005, 2012)

rm(democratic_exp)

#### Social Group Base Regime ####

table(wvsdata$regime_support_group_max) 


### delete missing macro-level cases 

wvsdata <- wvsdata %>%
  drop_na(e_migdppcln)

##################################
#Group and grand mean centering
##################################

# Author: Zoltan Fazekas
centFUN <- function(x) {
  x - mean(x, na.rm = TRUE)
}

# Median Cenetering for factor variables
centFUN_median <- function(x) {
  x - median(x, na.rm = TRUE)
}
#### Data manipulation ####

wvsdata <- wvsdata %>%
  mutate(democratic_expericene_all = ifelse(is.na(democratic_expericene_all), 0, democratic_expericene_all), 
         democratic_expericene = ifelse(is.na(democratic_expericene), 0, democratic_expericene), 
         regime_support_group_max = ifelse(is.na(regime_support_group_max), 0, regime_support_group_max))

# Use the function to grand-mean center the macro-level predictor
wvsdata$gini_mktCGM <- centFUN(wvsdata$gini_mkt)
wvsdata$gini_dispCGM <- centFUN(wvsdata$gini_disp)
wvsdata$v2csreprssCGM <- centFUN(wvsdata$v2csreprss)
wvsdata$e_migdppclnCGM <- centFUN(wvsdata$e_migdppcln)
wvsdata$regime_support_group_maxCGM <- centFUN(wvsdata$regime_support_group_max)
wvsdata$democratic_expericene_allCGM <- centFUN(wvsdata$democratic_expericene_all)
wvsdata$democratic_expericeneCGM <- centFUN(wvsdata$democratic_expericene)

summary(wvsdata$e_migdppclnCGM)
summary(wvsdata$v2csreprssCGM)
summary(wvsdata$gini_mktCGM)
summary(wvsdata$gini_dispCGM)
summary(wvsdata$regime_support_group_maxCGM)
summary(wvsdata$democratic_expericeneCGM)
summary(wvsdata$democratic_expericene_allCGM)

#Use the fuction to group-mean center the individual predictors
wvsdata <- wvsdata %>%
  group_by(S025) %>%
  mutate(ageCWC = centFUN(age))
summary(wvsdata$ageCWC)

wvsdata <- wvsdata %>%
  group_by(S025) %>%
  mutate(educationCWC = centFUN(as.numeric(education)),
         income_quintilesCWC = centFUN(as.numeric(income_quintiles)),
         income_decilesCWC = centFUN(as.numeric(income_deciles)))

summary(wvsdata$educationCWC)
summary(wvsdata$income_quintilesCWC)
summary(wvsdata$income_decilesCWC)

#Re-scaling predictor variables
wvsdata$ageCWC <- wvsdata$ageCWC/10


##################################
#Analysis: three-level models: Market Inequality
##################################

v1 <- glmer(civilsociety_b ~ 1 + (1 | S025) + (1 | S003) ,
            data = wvsdata, 
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v1)
icc(v1, adjusted = TRUE)
icc(v1)


v2 <- glmer(civilsociety_b ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_decilesCWC + (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v2)
anova(v1, v2)

#macrolevel gini
v3 <- glmer(civilsociety_b ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_decilesCWC + 
              gini_mktCGM + 
              (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v3)
anova(v2, v3)

#macrolevel gini plus addtional macro controls
v4 <- glmer(civilsociety_b ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_decilesCWC + 
              gini_mktCGM + v2csreprssCGM + e_migdppclnCGM + 
              democratic_expericeneCGM + regime_support_group_maxCGM +
              (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v4)
anova(v2, v4)

# model 4 plus random slopes
v5 <- glmer(civilsociety_b ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_decilesCWC + 
              gini_mktCGM + v2csreprssCGM + e_migdppclnCGM + 
              democratic_expericeneCGM + regime_support_group_maxCGM +
              (1 + income_decilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v5)
anova(v2, v5)

# model 5 plus interactions 
v6 <- glmer(civilsociety_b ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_decilesCWC +
              gini_mktCGM + v2csreprssCGM + e_migdppclnCGM + 
              democratic_expericeneCGM + regime_support_group_maxCGM +
              income_decilesCWC*gini_mktCGM + v2csreprssCGM*income_decilesCWC +
              (1 + income_decilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v6)
anova(v5, v6)


# model 6 but with interaction of freedom/fariness and income level and without gini
v7 <- glmer(civilsociety_b ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_decilesCWC +
              v2csreprssCGM + e_migdppclnCGM + 
              democratic_expericeneCGM + regime_support_group_maxCGM +
              income_decilesCWC*v2csreprssCGM +
              (1 + income_decilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v7)
anova(v6, v7)

##Interplot relationship between income and inequality

library(interplot)
plot_Interaction_civil_mkt <- interplot(v6, 
                                        var1 = "gini_mktCGM", 
                                        var2 = "income_decilesCWC", 
                                        ci = 0.95, 
                                        hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of inequality on civil society participation") +
  theme_bw() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "red",
             size = 0.5)
plot_Interaction_civil_mkt
ggsave("output/C3interactioninequality.pdf")


plot_Interaction_civil_repress <- interplot(v6, 
                                            var1 = "v2csreprssCGM", 
                                            var2 = "income_decilesCWC", 
                                            ci = 0.95, 
                                            hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of CSO repression on civil society participation") +
  theme_bw() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "red",
             size = 0.5)
plot_Interaction_civil_repress
ggsave("output/C3interactionrepressionv6.pdf")


plot_Interaction_civil_repress2 <- interplot(v7, 
                                             var1 = "v2csreprssCGM", 
                                             var2 = "income_decilesCWC", 
                                             ci = 0.95, 
                                             hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of CSO repression on civil society participation") +
  theme_bw() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "red",
             size = 0.5)
plot_Interaction_civil_repress2
ggsave("output/C3interactionrepressionv7.pdf")

######################################################
#Building Logg Odds ratios of Models

library(dotwhisker)
library(broom.mixed)

v1_df <- tidy(v1, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v2_df <- tidy(v2, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v3_df <- tidy(v3, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v4_df <- tidy(v4, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v5_df <- tidy(v5, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v6_df <- tidy(v6, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v7_df <- tidy(v7, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")

#Rescaling ordinal and continous predictors by tow standard deviations
v1_df <- v1_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 1") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_decilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_decilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_decilesCWC:v2csreprssCGM" = "Income * CSO repression"))

v2_df <- v2_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 2") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_decilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_decilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_decilesCWC:v2csreprssCGM" = "Income * CSO repression"))
v3_df <- v3_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 3") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_decilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_decilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_decilesCWC:v2csreprssCGM" = "Income * CSO repression"))
v4_df <- v4_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 4") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_decilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_decilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_decilesCWC:v2csreprssCGM" = "Income * CSO repression"))
v5_df <- v5_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 5") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_decilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_decilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_decilesCWC:v2csreprssCGM" = "Income * CSO repression"))

v6_df <- v6_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 6") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_decilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_decilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_decilesCWC:v2csreprssCGM" = "Income * CSO repression"))

v7_df <- v7_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 7") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_decilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_decilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_decilesCWC:v2csreprssCGM" = "Income * CSO repression"))

#Log Odds Plot Effects Plot

v_models <- rbind(v2_df, v3_df, v4_df, v5_df, v6_df, v7_df)

v_brackets <- list(c("Individual Predictors", "Male", "Income (centered)"), 
                   c("Country-year Predictors",  "Market Inequality (centered)", "Regime Support Group (centered)"),
                   c("Cross level Interactions", "Income * Inequality" , "Income * CSO repression"))

plot_robustness_civil <- {dwplot(v_models) +
    theme_bw() + xlab("Coefficient Estimate") + ylab("") +
    geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
    ggtitle("Predicting Civil Society Participation") +
    theme(plot.title = element_text(face="bold"),
          legend.position = c(0.8, 0.8),
          legend.justification = c(0, 0), 
          legend.background = element_rect(colour="grey80"),
          legend.title = element_blank())} %>% 
  add_brackets(v_brackets)

plot_robustness_civil
ggsave("output/C3dotwhisker.pdf", width = 22 , height = 25 , units = c("cm"))


#### texreg regression tables main models ####

library(texreg)
texreg(list(v1, v2, v3, v4, v5, v6, v7), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)","Male",
                             "Age(centered)","Education(centered)",
                             "Children",
                             "Unemployed",
                             "Income deciles (centered)",
                             "Inequality (centered)",
                             "CSO repression (centered",
                             "GDP pc log (centered)",
                             "Previous Dem. Experience (centered)", 
                             "Regime Support Group (centered)", 
                             "Income * Inequality",
                             "Income * CSO repression"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Civil Society Participation: Robustness Income Deciles",
       fontsize = "scriptsize",
       label = "Table:E1", 
       longtable = TRUE)

###############################
##Using Disposable Inequality##
###############################

d1 <- glmer(civilsociety_b ~ 1 + (1 | S025) + (1 | S003) ,
            data = wvsdata, 
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d1)
icc(d1, adjusted = TRUE)
icc(d1)


d2 <- glmer(civilsociety_b ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_decilesCWC + (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d2)
anova(d1, d2)

#macrolevel gini
d3 <- glmer(civilsociety_b ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_decilesCWC + 
              gini_dispCGM + 
              (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d3)
anova(d2, d3)

#macrolevel gini plus addtional macro controls
d4 <- glmer(civilsociety_b ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_decilesCWC + 
              gini_dispCGM + v2csreprssCGM + e_migdppclnCGM + 
              democratic_expericeneCGM + regime_support_group_maxCGM +
              (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d4)
anova(d2, d4)

# model 4 plus random slopes
d5 <- glmer(civilsociety_b ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_decilesCWC + 
              gini_dispCGM + v2csreprssCGM + e_migdppclnCGM + 
              democratic_expericeneCGM + regime_support_group_maxCGM +
              (1 + income_decilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d5)
anova(d2, d5)

# model 5 plus interactions 
d6 <- glmer(civilsociety_b ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_decilesCWC +
              gini_dispCGM + v2csreprssCGM + e_migdppclnCGM + 
              democratic_expericeneCGM + regime_support_group_maxCGM +
              income_decilesCWC*gini_dispCGM + v2csreprssCGM*income_decilesCWC +
              (1 + income_decilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d6)
anova(d5, d6)


# model 6 but with interaction of freedom/fariness and income level and without gini
d7 <- glmer(civilsociety_b ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_decilesCWC +
              v2csreprssCGM + e_migdppclnCGM + 
              democratic_expericeneCGM + regime_support_group_maxCGM +
              income_decilesCWC*v2csreprssCGM +
              (1 + income_decilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d7)
anova(d6, d7)

##Interplot relationship between income and inequality

library(interplot)
plot_Interaction_civil_disp <- interplot(d6, 
                                         var1 = "gini_dispCGM", 
                                         var2 = "income_decilesCWC", 
                                         ci = 0.95, 
                                         hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of inequality on civil society participation") +
  theme_bw() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "red",
             size = 0.5)
plot_Interaction_civil_disp
ggsave("output/C4interactioninequality.pdf")


plot_Interaction_civil_repress <- interplot(d6, 
                                            var1 = "v2csreprssCGM", 
                                            var2 = "income_decilesCWC", 
                                            ci = 0.95, 
                                            hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of CSO repression on civil society participation") +
  theme_bw() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "red",
             size = 0.5)
plot_Interaction_civil_repress
ggsave("output/C4interactionrepressionv6.pdf")


plot_Interaction_civil_repress2 <- interplot(d7, 
                                             var1 = "v2csreprssCGM", 
                                             var2 = "income_decilesCWC", 
                                             ci = 0.95, 
                                             hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of CSO repression on civil society participation") +
  theme_bw() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "red",
             size = 0.5)
plot_Interaction_civil_repress2
ggsave("output/C4interactionrepressionv7.pdf")

######################################################
#Building Logg Odds ratios of Models

library(dotwhisker)
library(broom.mixed)

d1_df <- tidy(d1, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
d2_df <- tidy(d2, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
d3_df <- tidy(d3, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
d4_df <- tidy(d4, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
d5_df <- tidy(d5, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
d6_df <- tidy(d6, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
d7_df <- tidy(d7, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")

#Rescaling ordinal and continous predictors by tow standard deviations
d1_df <- d1_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 1") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_decilesCWC = "Income (centered)",
                       gini_dispCGM = "Disp. Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_decilesCWC:gini_dispCGM" = "Income * Inequality",
                       "income_decilesCWC:v2csreprssCGM" = "Income * CSO repression"))

d2_df <- d2_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 2") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_decilesCWC = "Income (centered)",
                       gini_dispCGM = "Disp. Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_decilesCWC:gini_dispCGM" = "Income * Inequality",
                       "income_decilesCWC:v2csreprssCGM" = "Income * CSO repression"))

d3_df <- d3_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 3") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_decilesCWC = "Income (centered)",
                       gini_dispCGM = "Disp. Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_decilesCWC:gini_dispCGM" = "Income * Inequality",
                       "income_decilesCWC:v2csreprssCGM" = "Income * CSO repression"))

d4_df <- d4_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 4") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_decilesCWC = "Income (centered)",
                       gini_dispCGM = "Disp. Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_decilesCWC:gini_dispCGM" = "Income * Inequality",
                       "income_decilesCWC:v2csreprssCGM" = "Income * CSO repression"))

d5_df <- d5_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 5") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_decilesCWC = "Income (centered)",
                       gini_dispCGM = "Disp. Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_decilesCWC:gini_dispCGM" = "Income * Inequality",
                       "income_decilesCWC:v2csreprssCGM" = "Income * CSO repression"))

d6_df <- d6_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 6") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_decilesCWC = "Income (centered)",
                       gini_dispCGM = "Disp. Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_decilesCWC:gini_dispCGM" = "Income * Inequality",
                       "income_decilesCWC:v2csreprssCGM" = "Income * CSO repression"))

d7_df <- d7_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 7") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_decilesCWC = "Income (centered)",
                       gini_dispCGM = "Disp. Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_decilesCWC:gini_dispCGM" = "Income * Inequality",
                       "income_decilesCWC:v2csreprssCGM" = "Income * CSO repression"))

#Log Odds Plot Effects Plot

d_models <- rbind(d2_df, d3_df, d4_df, d5_df, d6_df, d7_df)

d_brackets <- list(c("Individual Predictors", "Male", "Income (centered)"), 
                   c("Country-year Predictors", "Disp. Inequality (centered)", "Regime Support Group (centered)"),
                   c("Cross level Interactions", "Income * Inequality" , "Income * CSO repression"))

plot_robustness_civil_disp <- {dwplot(d_models) +
    theme_bw() + xlab("Coefficient Estimate") + ylab("") +
    geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
    ggtitle("Predicting Civil Society Participation") +
    theme(plot.title = element_text(face="bold"),
          legend.position = c(0.8, 0.8),
          legend.justification = c(0, 0), 
          legend.background = element_rect(colour="grey80"),
          legend.title = element_blank())} %>% 
  add_brackets(d_brackets)

plot_robustness_civil_disp
ggsave("output/C4dotwhisker.pdf", width = 22 , height = 25 , units = c("cm"))


#### texreg regression tables main models ####

library(texreg)
texreg(list(d1, d2, d3, d4, d5, d6, d7), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)","Male",
                             "Age(centered)","Education(centered)",
                             "Children",
                             "Unemployed",
                             "Income deciles (centered)",
                             "Inequality (centered)",
                             "CSO repression (centered",
                             "GDP pc log (centered)",
                             "Previous Dem. Experience (centered)", 
                             "Regime Support Group (centered)", 
                             "Income * Inequality",
                             "Income * CSO repression"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Civil Society Participation: Robustess Income deciles",
       fontsize = "scriptsize",
       label = "Table:E2", 
       longtable = TRUE)

###############################################################################################
###############################################################################################
#### 2. All the Ginis Data 
#### Supplementary Appendix E3


# clear workspace
rm(list=ls())

# set working directory
setwd("M:/Promotion/Paper/WVS Paper Inequality and Voting in Autocracies/calculations_new")

# loading packages

library(countrycode)
library(tidyverse)
library(viridis)
library(haven)
library(lme4)
library(sjstats)
library(texreg)

#### Import Data ####
wvsdata <- readRDS("data/wvsdata.rds") 

#### Import all the Ginis measure by Vdem version 8.0 ####

vdem8 <- read_csv("data/vdem_8/V-Dem-CY+Others-v8.csv")

vdem8 <- vdem8 %>%
  dplyr::select(country_id, year, e_peginiwi)

#### Merge All the ginis with WVSdata

wvsdata <- wvsdata %>%
  left_join(vdem8, c("country_id", "year"))

rm(vdem8)

#### Reducing to authoritarian regimes ####

wvsdata <- wvsdata %>%
  filter(v2x_regime_amb<=4)

wvsdata <- wvsdata %>%
  mutate(country_year = stringr::str_c(country_name, year, sep=" "))

table(wvsdata$country_year, is.na(wvsdata$e_peginiwi))
#Replacing missings in Dependent Variables

wvsdata <- wvsdata %>%
  filter(A064>=0 | A066>=0 | A067>=0 | A068>=0 | A069>=0 | A070>=0 |
           A072>=0 | A073>=0 | A074>=0 |A075>=0 | A076>=0 | A077>=0
         | A078>=0 | A079>=0 | A099>=0 | A100>=0 | A101>=0 | A102>=0
         | A103>=0 | A104>=0 | A105>=0 | A106>=0 | A106B>=0 | A106C>=0 )%>%
  mutate(civilsociety_b = if_else(A064==1, 1,
                                  if_else(A066==1, 1,
                                          if_else(A067==1, 1, 
                                                  if_else(A068==1, 1,
                                                          if_else(A069==1, 1, 
                                                                  if_else(A070==1, 1,
                                                                          if_else(A072==1, 1, 
                                                                                  if_else(A073==1, 1,
                                                                                          if_else(A074==1, 1,
                                                                                                  if_else(A075==1, 1,
                                                                                                          if_else(A076==1, 1, 
                                                                                                                  if_else(A077==1, 1,
                                                                                                                          if_else(A078==1, 1,
                                                                                                                                  if_else(A079==1, 1, 
                                                                                                                                          if_else(A099>=1, 1,
                                                                                                                                                  if_else(A100>=1, 1, 
                                                                                                                                                          if_else(A101>=1, 1,
                                                                                                                                                                  if_else(A102>=1, 1,
                                                                                                                                                                          if_else(A103>=1, 1,
                                                                                                                                                                                  if_else(A104>=1, 1, 
                                                                                                                                                                                          if_else(A105>=1, 1,
                                                                                                                                                                                                  if_else(A106>=1, 1,
                                                                                                                                                                                                          if_else(A106B>=1, 1, 
                                                                                                                                                                                                                  if_else(A106C>=1, 1, 0)))))))))))))))))))))))))

table(wvsdata$civilsociety_b)

#Repalcing missings in Independent Variables

wvsdata <- wvsdata %>%
  filter(X001>=0)%>%
  mutate(sex= X001)

wvsdata <- wvsdata %>%
  filter(X003>=0)%>%
  mutate(age= X003)

wvsdata <- wvsdata %>%
  filter(X025>=0)%>%
  mutate(education= X025)

wvsdata <- wvsdata %>%
  filter(X028>=0)%>%
  mutate(employement_status= X028)

wvsdata <- wvsdata %>%
  filter(X047>=0)%>%
  mutate(income_deciles= X047)

wvsdata <- wvsdata %>%
  mutate(income_quintiles = ifelse(income_deciles==1, 1,
                                   ifelse(income_deciles==2, 1,
                                          ifelse(income_deciles==3, 2,
                                                 ifelse(income_deciles==4, 2,
                                                        ifelse(income_deciles==5, 3,
                                                               ifelse(income_deciles==6, 3,
                                                                      ifelse(income_deciles==7 , 4,
                                                                             ifelse(income_deciles==8, 4, 5)))))))))


wvsdata <- wvsdata %>%
  filter(X011>=0)%>%
  mutate(children= ifelse(X011>0, 1, 0))

wvsdata <- wvsdata %>%
  filter(X028>=0)%>%
  mutate(unemployed= ifelse(X028==7, 1, 0))

table(wvsdata$civilsociety_b)

# 35907 non member of civil society/ not belonged to civil society, 29471 member/ Belong to cso

###Data Manipulation of GINI measures, fill some missing value with last national observations

wvsdata <- wvsdata %>%
  drop_na(e_peginiwi)

table(wvsdata$civilsociety_b)
# 33529 non member of civil society/ not belonged to civil society, 26631 member/ Belong to cso

##################################
#Group and grand mean centering
##################################

# Author: Zoltan Fazekas
centFUN <- function(x) {
  x - mean(x, na.rm = TRUE)
}

# Median Centering for factor variables
centFUN_median <- function(x) {
  x - median(x, na.rm = TRUE)
}

#### Data manipulation ####

wvsdata <- wvsdata %>%
  mutate(democratic_expericene_all = ifelse(is.na(democratic_expericene_all), 0, democratic_expericene_all), 
         democratic_expericene = ifelse(is.na(democratic_expericene), 0, democratic_expericene), 
         regime_support_group_max = ifelse(is.na(regime_support_group_max), 0, regime_support_group_max))

# Use the function to grand-mean center the macro-level predictor
wvsdata$e_peginiwiCGM <- centFUN(wvsdata$e_peginiwi)
wvsdata$v2csreprssCGM <- centFUN(wvsdata$v2csreprss)
wvsdata$e_migdppclnCGM <- centFUN(wvsdata$e_migdppcln)
wvsdata$regime_support_group_maxCGM <- centFUN(wvsdata$regime_support_group_max)
wvsdata$democratic_expericene_allCGM <- centFUN(wvsdata$democratic_expericene_all)
wvsdata$democratic_expericeneCGM <- centFUN(wvsdata$democratic_expericene)


summary(wvsdata$e_migdppclnCGM)
summary(wvsdata$v2csreprssCGM)
summary(wvsdata$e_peginiwiCGM)
summary(wvsdata$regime_support_group_maxCGM)
summary(wvsdata$democratic_expericeneCGM)
summary(wvsdata$democratic_expericene_allCGM)


#Use the function to group-mean center the individual predictors
wvsdata <- wvsdata %>%
  group_by(S025) %>%
  mutate(ageCWC = centFUN(age))
summary(wvsdata$ageCWC)

wvsdata <- wvsdata %>%
  group_by(S025) %>%
  mutate(educationCWC = centFUN(as.numeric(education)),
         income_quintilesCWC = centFUN(as.numeric(income_quintiles)),
         income_decilesCWC = centFUN(as.numeric(income_deciles)))

summary(wvsdata$educationCWC)
summary(wvsdata$income_quintilesCWC)
summary(wvsdata$income_decilesCWC)

#Re-scaling predictor variables
wvsdata$ageCWC <- wvsdata$ageCWC/10

##################################
#Analysis: UNU WIDER GINI
##################################

v1 <- glmer(civilsociety_b ~ 1 + (1 | S025) + (1 | S003) ,
            data = wvsdata, 
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v1)
icc(v1, adjusted = TRUE)
icc(v1)


v2 <- glmer(civilsociety_b ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC + (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v2)
anova(v1, v2)

#macrolevel gini
v3 <- glmer(civilsociety_b ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC + 
              e_peginiwiCGM + 
              (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v3)
anova(v2, v3)

#macrolevel gini plus addtional macro controls
v4 <- glmer(civilsociety_b ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC + 
              e_peginiwiCGM + v2csreprssCGM + e_migdppclnCGM + 
              democratic_expericeneCGM + regime_support_group_maxCGM +
              (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v4)
anova(v2, v4)

# model 4 plus random slopes
v5 <- glmer(civilsociety_b ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC + 
              e_peginiwiCGM + v2csreprssCGM + e_migdppclnCGM + 
              democratic_expericeneCGM + regime_support_group_maxCGM +
              (1 + income_quintilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v5)
anova(v2, v5)

# model 5 plus interactions 
v6 <- glmer(civilsociety_b ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC +
              e_peginiwiCGM + v2csreprssCGM + e_migdppclnCGM + 
              democratic_expericeneCGM + regime_support_group_maxCGM +
              income_quintilesCWC*e_peginiwiCGM + v2csreprssCGM*income_quintilesCWC +
              (1 + income_quintilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v6)
anova(v5, v6)


# model 6 but with interaction of freedom/fariness and income level and without gini
v7 <- glmer(civilsociety_b ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC +
              v2csreprssCGM + e_migdppclnCGM + 
              democratic_expericeneCGM + regime_support_group_maxCGM +
              income_quintilesCWC*v2csreprssCGM +
              (1 + income_quintilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v7)
anova(v6, v7)

##Interplot relationship between income and inequality

library(interplot)
plot_Interaction_civil_mkt <- interplot(v6, 
                                        var1 = "e_peginiwiCGM", 
                                        var2 = "income_quintilesCWC", 
                                        ci = 0.95, 
                                        hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of inequality on civil society participation") +
  theme_bw() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "red",
             size = 0.5)
plot_Interaction_civil_mkt
ggsave("output/C5interactioninequality.pdf")


plot_Interaction_civil_repress <- interplot(v6, 
                                            var1 = "v2csreprssCGM", 
                                            var2 = "income_quintilesCWC", 
                                            ci = 0.95, 
                                            hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of CSO repression on civil society participation") +
  theme_bw() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "red",
             size = 0.5)
plot_Interaction_civil_repress
ggsave("output/C5interactionrepressionv6.pdf")


plot_Interaction_civil_repress2 <- interplot(v7, 
                                             var1 = "v2csreprssCGM", 
                                             var2 = "income_quintilesCWC", 
                                             ci = 0.95, 
                                             hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of CSO repression on civil society participation") +
  theme_bw() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "red",
             size = 0.5)
plot_Interaction_civil_repress2
ggsave("output/C5interactionrepressionv7.pdf")

######################################################
#Building Logg Odds ratios of Models

library(dotwhisker)
library(broom.mixed)

v1_df <- tidy(v1, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v2_df <- tidy(v2, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v3_df <- tidy(v3, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v4_df <- tidy(v4, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v5_df <- tidy(v5, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v6_df <- tidy(v6, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v7_df <- tidy(v7, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")

#Rescaling ordinal and continous predictors by tow standard deviations
v1_df <- v1_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 1") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       e_peginiwiCGM = "Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:e_peginiwiCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2csreprssCGM" = "Income * CSO repression"))

v2_df <- v2_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 2") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       e_peginiwiCGM = "Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:e_peginiwiCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2csreprssCGM" = "Income * CSO repression"))

v3_df <- v3_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 3") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       e_peginiwiCGM = "Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:e_peginiwiCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2csreprssCGM" = "Income * CSO repression"))

v4_df <- v4_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 4") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       e_peginiwiCGM = "Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:e_peginiwiCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2csreprssCGM" = "Income * CSO repression"))

v5_df <- v5_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 5") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       e_peginiwiCGM = "Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:e_peginiwiCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2csreprssCGM" = "Income * CSO repression"))


v6_df <- v6_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 6") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       e_peginiwiCGM = "Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:e_peginiwiCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2csreprssCGM" = "Income * CSO repression"))


v7_df <- v7_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 7") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       e_peginiwiCGM = "Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:e_peginiwiCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2csreprssCGM" = "Income * CSO repression"))

#Log Odds Plot Effects Plot

v_models <- rbind(v2_df, v3_df, v4_df, v5_df, v6_df, v7_df)

v_brackets <- list(c("Individual Predictors", "Male", "Income (centered)"), 
                   c("Country-year Predictors",  "Inequality (centered)", "Regime Support Group (centered)"),
                   c("Cross level Interactions", "Income * Inequality" , "Income * CSO repression"))

plot_robust_civil_unu <- {dwplot(v_models) +
    theme_bw() + xlab("Coefficient Estimate") + ylab("") +
    geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
    ggtitle("Predicting Civil Society Participation") +
    theme(plot.title = element_text(face="bold"),
          legend.position = c(0.8, 0.8),
          legend.justification = c(0, 0), 
          legend.background = element_rect(colour="grey80"),
          legend.title = element_blank())} %>% 
  add_brackets(v_brackets)

plot_robust_civil_unu
ggsave("output/C5dotwhisker.pdf", width = 22 , height = 25 , units = c("cm"))


#### texreg regression tables main models ####

library(texreg)
texreg(list(v1, v2, v3, v4, v5, v6, v7), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)","Male",
                             "Age(centered)","Education(centered)",
                             "Children",
                             "Unemployed",
                             "Income quintiles (centered)",
                             "Inequality (centered)",
                             "CSO repression (centered)",
                             "GDP pc log (centered)",
                             "Previous Dem. Experience (centered)", 
                             "Regime Support Group (centered)", 
                             "Income * Inequality",
                             "Income * CSO repression"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Civil Society Particpation. Robustness UNU WIDER GINI",
       fontsize = "scriptsize",
       label = "Table:E4", 
       longtable = TRUE)

###############################################################################################
###############################################################################################
#### 3. Civil Society variable without party-near organizations
#### Supplementary Appendix E4


#### Import Data ####
wvsdata <- readRDS("data/wvsdata.rds") 

#### Reducing to authoritarian regimes ####

wvsdata <- wvsdata %>%
  filter(v2x_regime_amb<=4)

wvsdata <- wvsdata %>%
  mutate(country_year = stringr::str_c(country_name, year, sep=" "))

table(wvsdata$country_year, is.na(wvsdata$gini_mkt))

#### Building variables ####

# Dependent variable without party organizsations 
wvsdata <- wvsdata %>%
  filter(A064>=0 | A066>=0 | A067>=0 | A069>=0 | A070>=0 |
           A072>=0 | A073>=0 | A074>=0 | A075>=0 | A076>=0 | A077>=0
         | A078>=0 | A079>=0 | A099>=0 | A100>=0 | A101>=0 
         | A103>=0 | A104>=0 | A105>=0 | A106>=0 | A106B>=0 | A106C>=0 ) %>%
  mutate(civilsociety_p = if_else(A064==1, 1,
                                  if_else(A066==1, 1,
                                          if_else(A067==1, 1, 
                                                  if_else(A069==1, 1, 
                                                          if_else(A070==1, 1,
                                                                  if_else(A072==1, 1, 
                                                                          if_else(A073==1, 1,
                                                                                  if_else(A074==1, 1,
                                                                                          if_else(A075==1, 1,
                                                                                                  if_else(A076==1, 1, 
                                                                                                          if_else(A077==1, 1,
                                                                                                                  if_else(A078==1, 1,
                                                                                                                          if_else(A079==1, 1, 
                                                                                                                                  if_else(A099>=1, 1,
                                                                                                                                          if_else(A100>=1, 1, 
                                                                                                                                                  if_else(A101>=1, 1,
                                                                                                                                                          if_else(A103>=1, 1,
                                                                                                                                                                  if_else(A104>=1, 1, 
                                                                                                                                                                          if_else(A105>=1, 1,
                                                                                                                                                                                  if_else(A106>=1, 1,
                                                                                                                                                                                          if_else(A106B>=1, 1, 
                                                                                                                                                                                                  if_else(A106C>=1, 1, 0)))))))))))))))))))))))

table(wvsdata$civilsociety_p)

#Repalcing missings in Independent Variables

wvsdata <- wvsdata %>%
  filter(X001>=0)%>%
  mutate(sex= X001)

wvsdata <- wvsdata %>%
  filter(X003>=0)%>%
  mutate(age= X003)

wvsdata <- wvsdata %>%
  filter(X025>=0)%>%
  mutate(education= X025)

wvsdata <- wvsdata %>%
  filter(X028>=0)%>%
  mutate(employement_status= X028)

wvsdata <- wvsdata %>%
  filter(X047>=0)%>%
  mutate(income_deciles= X047)

wvsdata <- wvsdata %>%
  mutate(income_quintiles = ifelse(income_deciles==1, 1,
                                   ifelse(income_deciles==2, 1,
                                          ifelse(income_deciles==3, 2,
                                                 ifelse(income_deciles==4, 2,
                                                        ifelse(income_deciles==5, 3,
                                                               ifelse(income_deciles==6, 3,
                                                                      ifelse(income_deciles==7 , 4,
                                                                             ifelse(income_deciles==8, 4, 5)))))))))


wvsdata <- wvsdata %>%
  filter(X011>=0)%>%
  mutate(children= ifelse(X011>0, 1, 0))

wvsdata <- wvsdata %>%
  filter(X028>=0)%>%
  mutate(unemployed= ifelse(X028==7, 1, 0))

table(wvsdata$civilsociety_p)

# 37703 non member of civil society/ not belonged to civil society, 27632 member/ Belong to cso

###Data Manipulation of GINI measures, fill some missing value with last national observations

missing <- wvsdata %>% 
  group_by(country_year) %>% 
  summarise(sumNA = sum(is.na(gini_mkt)))
missing

wvsdata <- wvsdata %>%
  mutate(gini_mkt=ifelse(country_year=="Algeria 2014", 36.7, gini_mkt)) %>% # Algeria 2011
  mutate(gini_disp=ifelse(country_year=="Algeria 2014", 35, gini_disp)) %>%
  mutate(gini_mkt=ifelse(country_year=="Azerbaijan 2011", 38.3, gini_mkt)) %>% # Azerbaijan 2011
  mutate(gini_disp=ifelse(country_year=="Azerbaijan 2011", 30.3, gini_disp)) %>%
  mutate(gini_disp=ifelse(country_year=="Lebanon 2013", 38.5, gini_mkt)) %>% # Lebanon 2012
  mutate(gini_disp=ifelse(country_year=="Lebanon 2013", 36.5, gini_disp))

wvsdata <- wvsdata %>%
  mutate(gini_mkt=ifelse(country_year=="Haiti 2016", 57.6, gini_mkt)) %>%
  mutate(gini_disp=ifelse(country_year=="Haiti 2016", 53.7, gini_disp)) %>%
  mutate(gini_mkt=ifelse(country_year=="Nigeria 2012", 46.8, gini_mkt)) %>%
  mutate(gini_disp=ifelse(country_year=="Nigeria 2012", 44.1, gini_disp)) 

### delete Kuwait 2014, Uzbekistan 2011 and Libya 2014, no GINI information in a range of maximum 4 years before or after the country-year observations

wvsdata <- wvsdata %>%
  drop_na(gini_mkt)
rm(missing)

table(wvsdata$civilsociety_b)
# 36955 non member of civil society/ not belonged to civil society, 30677 member/ Belong to cso


##################################
#Group and grand mean centering
##################################

# Author: Zoltan Fazekas
centFUN <- function(x) {
  x - mean(x, na.rm = TRUE)
}

# Median Cenetering for factor variables
centFUN_median <- function(x) {
  x - median(x, na.rm = TRUE)
}

#### Data manipulation ####

wvsdata <- wvsdata %>%
  mutate(democratic_expericene_all = ifelse(is.na(democratic_expericene_all), 0, democratic_expericene_all), 
         democratic_expericene = ifelse(is.na(democratic_expericene), 0, democratic_expericene), 
         regime_support_group_max = ifelse(is.na(regime_support_group_max), 0, regime_support_group_max))

wvsdata <- wvsdata %>%
  drop_na(e_migdppcln)

# Use the function to grand-mean center the macro-level predictor
wvsdata$gini_mktCGM <- centFUN(wvsdata$gini_mkt)
wvsdata$gini_dispCGM <- centFUN(wvsdata$gini_disp)
wvsdata$v2csreprssCGM <- centFUN(wvsdata$v2csreprss)
wvsdata$e_migdppclnCGM <- centFUN(wvsdata$e_migdppcln)
wvsdata$regime_support_group_maxCGM <- centFUN(wvsdata$regime_support_group_max)
wvsdata$democratic_expericene_allCGM <- centFUN(wvsdata$democratic_expericene_all)
wvsdata$democratic_expericeneCGM <- centFUN(wvsdata$democratic_expericene)

summary(wvsdata$e_migdppclnCGM)
summary(wvsdata$v2csreprssCGM)
summary(wvsdata$gini_mktCGM)
summary(wvsdata$gini_dispCGM)
summary(wvsdata$regime_support_group_maxCGM)
summary(wvsdata$democratic_expericeneCGM)
summary(wvsdata$democratic_expericene_allCGM)


#Use the function to group-mean center the individual predictors
wvsdata <- wvsdata %>%
  group_by(S025) %>%
  mutate(ageCWC = centFUN(age))
summary(wvsdata$ageCWC)

wvsdata <- wvsdata %>%
  group_by(S025) %>%
  mutate(educationCWC = centFUN(as.numeric(education)),
         income_quintilesCWC = centFUN(as.numeric(income_quintiles)),
         income_decilesCWC = centFUN(as.numeric(income_deciles)))

summary(wvsdata$educationCWC)
summary(wvsdata$income_quintilesCWC)
summary(wvsdata$income_decilesCWC)

#Rescaling predicator variables
wvsdata$ageCWC <- wvsdata$ageCWC/10


##################################
#Analysis: three-level models: Market Inequality
##################################

v1 <- glmer(civilsociety_p ~ 1 + (1 | S025) + (1 | S003) ,
            data = wvsdata, 
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v1)
icc(v1, adjusted = TRUE)
icc(v1)


v2 <- glmer(civilsociety_p ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC + (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v2)
anova(v1, v2)

#macrolevel gini
v3 <- glmer(civilsociety_p ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC + 
              gini_mktCGM + 
              (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v3)
anova(v2, v3)

#macro-level gini plus additional macro controls
v4 <- glmer(civilsociety_p ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC + 
              gini_mktCGM + v2csreprssCGM + e_migdppclnCGM + 
              democratic_expericeneCGM + regime_support_group_maxCGM +
              (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v4)
anova(v2, v4)

# model 4 plus random slopes
v5 <- glmer(civilsociety_p ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC + 
              gini_mktCGM + v2csreprssCGM + e_migdppclnCGM + 
              democratic_expericeneCGM + regime_support_group_maxCGM +
              (1 + income_quintilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v5)
anova(v2, v5)

# model 5 plus interactions 
v6 <- glmer(civilsociety_p ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC +
              gini_mktCGM + v2csreprssCGM + e_migdppclnCGM + 
              democratic_expericeneCGM + regime_support_group_maxCGM +
              income_quintilesCWC*gini_mktCGM + v2csreprssCGM*income_quintilesCWC +
              (1 + income_quintilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v6)
anova(v5, v6)


# model 6 but with interaction of freedom/fairness and income level and without gini
v7 <- glmer(civilsociety_p ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC +
              v2csreprssCGM + e_migdppclnCGM + 
              democratic_expericeneCGM + regime_support_group_maxCGM +
              income_quintilesCWC*v2csreprssCGM +
              (1 + income_quintilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(v7)
anova(v6, v7)

##Interplot relationship between income and inequality

library(interplot)
plot_Interaction_civil_mkt <- interplot(v6, 
                                        var1 = "gini_mktCGM", 
                                        var2 = "income_quintilesCWC", 
                                        ci = 0.95, 
                                        hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of inequality on civil society participation") +
  theme_bw() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "red",
             size = 0.5)
plot_Interaction_civil_mkt
ggsave("output/C6interactioninequality.pdf")


plot_Interaction_civil_repress <- interplot(v6, 
                                            var1 = "v2csreprssCGM", 
                                            var2 = "income_quintilesCWC", 
                                            ci = 0.95, 
                                            hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of CSO repression on civil society participation") +
  theme_bw() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "red",
             size = 0.5)
plot_Interaction_civil_repress
ggsave("output/C6interactionrepressionv6.pdf")


plot_Interaction_civil_repress2 <- interplot(v7, 
                                             var1 = "v2csreprssCGM", 
                                             var2 = "income_quintilesCWC", 
                                             ci = 0.95, 
                                             hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of CSO repression on civil society participation") +
  theme_bw() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "red",
             size = 0.5)
plot_Interaction_civil_repress2
ggsave("output/C6interactionrepressionv7.pdf")

######################################################
#Building Logg Odds ratios of Models

library(dotwhisker)
library(broom.mixed)

v1_df <- tidy(v1, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v2_df <- tidy(v2, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v3_df <- tidy(v3, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v4_df <- tidy(v4, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v5_df <- tidy(v5, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v6_df <- tidy(v6, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
v7_df <- tidy(v7, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")

#Rescaling ordinal and continous predictors by tow standard deviations
v1_df <- v1_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 1") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2csreprssCGM" = "Income * CSO repression"))

v2_df <- v2_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 2") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2csreprssCGM" = "Income * CSO repression"))
v3_df <- v3_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 3") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2csreprssCGM" = "Income * CSO repression"))
v4_df <- v4_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 4") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2csreprssCGM" = "Income * CSO repression"))
v5_df <- v5_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 5") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2csreprssCGM" = "Income * CSO repression"))

v6_df <- v6_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 6") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2csreprssCGM" = "Income * CSO repression"))

v7_df <- v7_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 7") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_mktCGM = "Market Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_mktCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2csreprssCGM" = "Income * CSO repression"))
#Log Odds Plot Effects Plot

v_models <- rbind(v2_df, v3_df, v4_df, v5_df, v6_df, v7_df)

v_brackets <- list(c("Individual Predictors", "Male", "Income (centered)"), 
                   c("Country-year Predictors",  "Market Inequality (centered)", "Regime Support Group (centered)"),
                   c("Cross level Interactions", "Income * Inequality" , "Income * CSO repression"))

plot_robustness_civil <- {dwplot(v_models) +
    theme_bw() + xlab("Coefficient Estimate") + ylab("") +
    geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
    ggtitle("Predicting Civil Society Participation") +
    theme(plot.title = element_text(face="bold"),
          legend.position = c(0.8, 0.8),
          legend.justification = c(0, 0), 
          legend.background = element_rect(colour="grey80"),
          legend.title = element_blank())} %>% 
  add_brackets(v_brackets)

plot_robustness_civil
ggsave("output/C6dotwhisker.pdf", width = 22 , height = 25 , units = c("cm"))


#### texreg regression tables main models ####

library(texreg)
texreg(list(v1, v2, v3, v4, v5, v6, v7), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)","Male",
                             "Age(centered)","Education(centered)",
                             "Children",
                             "Unemployed",
                             "Income quintiles (centered)",
                             "Inequality (centered)",
                             "CSO repression (centered",
                             "GDP pc log (centered)",
                             "Previous Dem. Experience (centered)", 
                             "Regime Support Group (centered)", 
                             "Income * Inequality",
                             "Income * CSO repression"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Main Models Voting",
       fontsize = "scriptsize",
       label = "Table:E5", 
       longtable = TRUE)

###############################
##Using Disposable Inequality##
################################

#### Supplementary Appendix E5 

d1 <- glmer(civilsociety_p ~ 1 + (1 | S025) + (1 | S003) ,
            data = wvsdata, 
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d1)
icc(d1, adjusted = TRUE)
icc(d1)


d2 <- glmer(civilsociety_p ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC + (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d2)
anova(d1, d2)

#macrolevel gini
d3 <- glmer(civilsociety_p ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC + 
              gini_dispCGM + 
              (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d3)
anova(d2, d3)

#macrolevel gini plus addtional macro controls
d4 <- glmer(civilsociety_p ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC + 
              gini_dispCGM + v2csreprssCGM + e_migdppclnCGM + 
              democratic_expericeneCGM + regime_support_group_maxCGM +
              (1 | S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d4)
anova(d2, d4)

# model 4 plus random slopes
d5 <- glmer(civilsociety_p ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC + 
              gini_dispCGM + v2csreprssCGM + e_migdppclnCGM + 
              democratic_expericeneCGM + regime_support_group_maxCGM +
              (1 + income_quintilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d5)
anova(d2, d5)

# model 5 plus interactions 
d6 <- glmer(civilsociety_p ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC +
              gini_dispCGM + v2csreprssCGM + e_migdppclnCGM + 
              democratic_expericeneCGM + regime_support_group_maxCGM +
              income_quintilesCWC*gini_dispCGM + v2csreprssCGM*income_quintilesCWC +
              (1 + income_quintilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d6)
anova(d5, d6)


# model 6 but with interaction of freedom/fariness and income level and without gini
d7 <- glmer(civilsociety_p ~ sex + ageCWC+ educationCWC + children + 
              unemployed + income_quintilesCWC +
              v2csreprssCGM + e_migdppclnCGM + 
              democratic_expericeneCGM + regime_support_group_maxCGM +
              income_quintilesCWC*v2csreprssCGM +
              (1 + income_quintilesCWC| S025) + (1 | S003),
            data = wvsdata,
            na.action = na.omit,
            family = binomial(link = "logit"),
            control=glmerControl(optCtrl=list(maxfun=6000)))
summary(d7)
anova(d6, d7)

##Interplot relationship between income and inequality

library(interplot)
plot_Interaction_civil_disp <- interplot(d6, 
                                         var1 = "gini_dispCGM", 
                                         var2 = "income_quintilesCWC", 
                                         ci = 0.95, 
                                         hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of inequality on civil society participation") +
  theme_bw() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "red",
             size = 0.5)
plot_Interaction_civil_disp
ggsave("output/C7interactioninequality.pdf")


plot_Interaction_civil_repress <- interplot(d6, 
                                            var1 = "v2csreprssCGM", 
                                            var2 = "income_quintilesCWC", 
                                            ci = 0.95, 
                                            hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of CSO repression on civil society participation") +
  theme_bw() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "red",
             size = 0.5)
plot_Interaction_civil_repress
ggsave("output/C7interactionrepressionv6.pdf")


plot_Interaction_civil_repress2 <- interplot(d7, 
                                             var1 = "v2csreprssCGM", 
                                             var2 = "income_quintilesCWC", 
                                             ci = 0.95, 
                                             hist = TRUE) + 
  xlab("Individual Income") +
  ylab("Effect of CSO repression on civil society participation") +
  theme_bw() +
  geom_hline(yintercept = 0, 
             linetype = "dashed",
             color = "red",
             size = 0.5)
plot_Interaction_civil_repress2
ggsave("output/C7interactionrepressionv7.pdf")

######################################################
#Building Logg Odds ratios of Models

library(dotwhisker)
library(broom.mixed)

d1_df <- tidy(d1, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
d2_df <- tidy(d2, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
d3_df <- tidy(d3, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
d4_df <- tidy(d4, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
d5_df <- tidy(d5, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
d6_df <- tidy(d6, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
d7_df <- tidy(d7, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")

#Rescaling ordinal and continous predictors by tow standard deviations
d1_df <- d1_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 1") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_dispCGM = "Disp. Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_dispCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2csreprssCGM" = "Income * CSO repression"))

d2_df <- d2_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 2") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_dispCGM = "Disp. Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_dispCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2csreprssCGM" = "Income * CSO repression"))

d3_df <- d3_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 3") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_dispCGM = "Disp. Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_dispCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2csreprssCGM" = "Income * CSO repression"))

d4_df <- d4_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 4") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_dispCGM = "Disp. Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_dispCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2csreprssCGM" = "Income * CSO repression"))

d5_df <- d5_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 5") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_dispCGM = "Disp. Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_dispCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2csreprssCGM" = "Income * CSO repression"))

d6_df <- d6_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 6") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_dispCGM = "Disp. Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_dispCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2csreprssCGM" = "Income * CSO repression"))

d7_df <- d7_df %>%
  by_2sd(wvsdata) %>%
  arrange(term) %>%
  mutate(model="Model 7") %>%
  relabel_predictors(c(sex = "Male",                       
                       ageCWC = "Age (centered)", 
                       educationCWC = "Education (centered)", 
                       children = "Children", 
                       unemployed = "Unemployed", 
                       income_quintilesCWC = "Income (centered)",
                       gini_dispCGM = "Disp. Inequality (centered)",
                       v2csreprssCGM = "CSO repression (centered)",
                       e_migdppclnCGM = "GDP pc log (centered)",
                       democratic_expericeneCGM = "Previos Dem. Experience (centered)", 
                       regime_support_group_maxCGM = "Regime Support Group (centered)",
                       "income_quintilesCWC:gini_dispCGM" = "Income * Inequality",
                       "income_quintilesCWC:v2csreprssCGM" = "Income * CSO repression"))

#Log Odds Plot Effects Plot

d_models <- rbind(d2_df, d3_df, d4_df, d5_df, d6_df, d7_df)

d_brackets <- list(c("Individual Predictors", "Male", "Income (centered)"), 
                   c("Country-year Predictors", "Disp. Inequality (centered)", "Regime Support Group (centered)"),
                   c("Cross level Interactions", "Income * Inequality" , "Income * CSO repression"))

plot_robust_civil_disp <- {dwplot(d_models) +
    theme_bw() + xlab("Coefficient Estimate") + ylab("") +
    geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
    ggtitle("Predicting Civil Society Participation") +
    theme(plot.title = element_text(face="bold"),
          legend.position = c(0.8, 0.8),
          legend.justification = c(0, 0), 
          legend.background = element_rect(colour="grey80"),
          legend.title = element_blank())} %>% 
  add_brackets(d_brackets)

plot_robust_civil_disp
ggsave("output/C7dotwhisker.pdf", width = 22 , height = 25 , units = c("cm"))


#### texreg regression tables main models ####

library(texreg)
texreg(list(d1, d2, d3, d4, d5, d6, d7), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)","Male",
                             "Age(centered)","Education(centered)",
                             "Children",
                             "Unemployed",
                             "Income quintiles (centered)",
                             "Inequality (centered)",
                             "CSO repression (centered",
                             "GDP pc log (centered)",
                             "Previous Dem. Experience (centered)", 
                             "Regime Support Group (centered)", 
                             "Income * Inequality",
                             "Income * CSO repression"),
       booktabs = TRUE,
       use.packages = FALSE,
       fontsize = "scriptsize",
       label = "Table:E5", 
       longtable = TRUE)
